Visualization of stacking faults in fee crystals in plastic deformations 
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Using molecular dynamics simulation, we investigate the dynamics of stacking faults in fee crystals 
in uniaxial stretching in a Lennard- Jones binary mixture composed of 4096 particles in three dimen- 
sions. We visualize stacking faults using a disorder variable Dj{t) for each particle j constructed 
from local bond order parameters based on spherical harmonics (Steinhardt order parameters). Also 
introducing a method of bond breakage, we examine how stacking faults are formed and removed 
by collective particle motions. These processes are relevant in plasticity of fee crystals. 

PACS numbers: 83.10.Rs, 61.72.Nn, 62.20.F-, 61.66.-f 



I. INTRODUCTION 

In crystal plasticity, dislocations play a major rol^SEl, 
They tend to appear around point defects, preexisting 
dislocations, grain boundaries, and precipitates. They 
then grow to form slip planes, which extend over large 
spatial regions under the strong influence of the crystal 
structure^Ki^, polycrystal, birth and growth of dislo- 
cations are strongly influenced by grain boundaries, so 
the grain size is a key factor of plasticity. With decreas- 
ing the grain size, sliding motions of the particles at grain 
boundaries become increasingly important!^, as h as been 
revealed by molecular dynamics simulation j'' l ^^l t^. In 
crystal and polycrystal, plastic events have been observed 
as bursts or avalanches spanning wide ranges of space and 
time scales. We may mention experiments of acoustic 
emissiorP and transmission electron microscopj^^ . 

In particular, fee crystals may be regarded as a se- 
quence of closed-packed {111} layerd^ symbolically rep- 
resented by ABC ABC ABC , while hep crystals are rep- 
resented by ABABAB. If the free energy of the fee struc- 
ture is slightly lower than that of the hep structure, the 
regular fee sequence is often violated by stacking faults 
or twin faults. For example, an intrinsic stacking fault 
is represented by ABC AC ABC , where B is missing be- 
tween the middle two-particle layers AC. This planar 
defect is produced by collective slip motions B ^ C, 
C ^ A, and A ^ B {B ^ A, C ^ B, and A ^ C) on 
the right (left) side of the defect plane. If its size is finite 
within the crystal, a partial dislocation can be found at 
its border with a Burgers vector of type (a/6) (121), where 
a is the lattice constant. As is well-knowrP^, this stacking 
fault can disappear with superposition of another partial 
displacement of type (a/6)(112), where a full dislocation 
of type (a/2) (Oil) remains at its border. Plastic defor- 
mations can be achieved by collective stacking reorgani- 
zation for metals with small grain sizes (nanocrystals)Pl^. 
Stability of partials depend on various factors such as the 
grain size, the applied stain rate, and the potential en- 
ergies for the partial displacement. Full dislocation slips 
have rather been observed in metals with large grain sizes 
and with relatively large stacking fault energies. 

In colloids, direct observations of particle configura- 
tions are particularly informativ^^SHSOl ^-^^ colloid 



interaction is short-ranged and hard-core-like, there is 
almost no free energy difference between the fee and 
hep stucturefl^, so fee and hep stacking layers often 
appear randomly as random hexagonal closed-packed 
(rhcp) states^SHSll, A pure fee crystal can be realized 
by sedimentation of colloid particles onto a patterned 
[100] substrat^^. It is known that the earth gravity and 
small mechanical perturbations strongly affect the struc- 
ture of colloid crystal^. In oscillatory shear, the spatial 
distribution of stacking faults was measurecP^. Disloca- 
tions and s lip pl anes have been observed in indentation 
experiments^SEni. 

In three dimensions, disordered crystal states are very 
complex under the strong infiuence of the crystal struc- 
ture. To analyze them, Steinhardt et alW^ introduced 
tensorial bond-order parameters based on the spherical 
harmonics. With their method in analysis, the crystal- 
lization process has been simulatecP^ltSZI Stacking faults 
have been observed in growing crystal domains^SiSSl Jt 
has also b een us ed to examine the structural heterogene- 
ity in glas d^^ l '^^ . In this paper, we will present results of 
molecular dynamics simulation on the fundamental dy- 
namical processes of stacking faults in uniaxial stretch- 
ing. In this problem, it is highly desirable to develop 
the method of visualizing these collective and complex 
phenomena. With this purpose, we will introduce a dis- 
order variable Dj constructed from the Steinhardt order 
parameterf{2SI, which represents the deviation of hexago- 
nal crystalline order for each particle j. In our simulation, 
we will domonstrate that appearance and disappearnce 
of intrinsic stacking faults give rise to crystal plasticity. 

The organization of this paper is as follows. In Sec. II, 
our simulation method will be explained, where a method 
of bond breakage and the disorder variable Dj [t) will be 
introduced. In Sec. Ill, we will visualize stacking faults 
on the basis of our simulation to understand the large- 
scale configuration changes in plastic flow. We will exam- 
ine rapid time-evolution of the potential energies and the 
displacements of the particles close to a stacking fault at 
plastic events. We shall see that the edge of an expand- 
ing stacking fault propagates with a velocity close to the 
transverse sound velocity, while the particle velocities at 
the slip plane are much slower. We will also display the 
averages of various physical quantities over all the parti- 
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cles, including the stress-strain curve, on long time scales. 
They exhibit intermittent changes upon plastic events. 



II. BACKGROUND OF SIMULATION 

A. Model and simulation method 

We treat three-dimensional binary mixtures composed 
of two atomic species 1 and 2. The total particle number 
is TV = iVi + iVa = 4096. The composition of the large 
particles is written as c = N2/N. The particles interact 
via truncated Lennard- Jones (LJ) potentials, 



r 



\ r 



Q/3, 



(1) 



which are characterized by the energy e and the interac- 
tion lengths ctq/s = {da +cri3)/'2' {oi,P = 1,2). Here cti 
and (T2 are the soft-core diameters of the two species. For 
f > fcut = 2.25f72, we set Vap = and the constant Cap 
ensures the continuity of Vap at r = r^ut- The mass ra- 
tio is fixed at mi/m.2 — {(^1/(^2)^ ■ Initially, the particles 
were in a cubic box, whose length Lq was determined by 
Ll = Nial + N2crl. For c = 1, we set Lq = 15.85(72. The 
time step of integration is 0.005ro with 



To = 0-2 
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After preparation of the initial particle configurations 
at r = 0.015e/A:B, we applied uniaxial strain along the z 
axis for i > 0. To this end, we moved the position of the 
top plate of the cell z = L(t) as 



L(t) = Lo(l + £t), 
where e is the elongation rate chosen to be 
e = 3.0 X lO^Vn"^ 



(3) 
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All the side walls parallel to the z axis were deformable 
free boundariespSl. This free-surface condition was possi- 
ble due to the attractive parts of the LJ potentials. As 
in our two-dimensional simulatiorP^I^ we divided the cell 
into a bottom layer in the region < z < dp = 0.0625io, 
a top layer in the region Lit) — do < z < L{t), and 
an interior region in the region dp < ^ < L{t) — d^. 
Each layer contained pinning centers fixed at positions 
Rj (j = 1, • • • , A^f,). At each pinning center, a particle at 
position Tj is bound by a spring potential. 



u,{r,-R,)=K,\r,-R,\'l2, 



(5) 



where we set Kg — 2Qea2^ . In this paper, the num- 
ber of the bound particles in each layer is chosen to be 
N}) = 256. Then the number of the unbound particles 
in the interior is A^ub = N — 2Nt, — 3584. The positions 
of the pinning centers were independent of time in the 
bottom layer, but depended on time in the top layer as 



Rj{t) = Rj{0) + LitEz, where is the unit vector in 
the z axis. The bound particles in the layers interacted 
with the other unbound and bound particles via the LJ 
potentials in Eq.(l). The total potential energy consists 
of two parts as 



I Va(i{\rj-rk\)+ Uj{\rj-Rj\), (6) 
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where the summation in the first term is over all the 
particle pairs and that in the second term is over the 
2Ni, bound particles. 

The Nuh unbound particles were governed by the usual 
Newton equations. Their positions rj obeyed 
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where rj — d^rj/dt^. In our simulation, the unbound 
particles in touch with the bound particles turned out 
to be clamped to the boundary layers without slip and 
detachment during stretching. They also did not pen- 
etrated into the layers deeper than 172 ■ On th e other 
hand, we attached a Nose-Hoover thermostalP^IlIl to each 
boundary layer. That is, we introduced two thermostat 
variables Cbot(^) and Ctop(i)- Let B represent the top or 
bottom layer. The bound particles j E B obeyed 
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where rj = drj/dt, and t)g is the boundary velocity equal 
to LeEz at the top and to at the bottom. The thermo- 
stat variables Ce were governed by 
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where tnh is the thermostat characteristic time. In this 
paper, we used a very short time tnh — 0.072to. Then 
the top and bottom layers played the role of efficient ther- 
mostats subtracting the extra energy released at plastic 
deformations in the bulk. As a result, the local temper- 
ature (the particle kinetic energy averaged in a narrow 
region) became nearly homogeneous except for short du- 
ration periods of plastic events in the bulk region. 

We explain our preparation method of the initial par- 
ticle configurations^^' We followed four steps in the 
time region t < by solving Eqs.(6) and (7) with vb = 0. 
In the first three steps, we imposed the periodic bound- 
ary condition in all the directions with N = 4096. (i) 
First, we created random particle configurations in a 
liquid state at T = l.Tbe/ks in a time interval longer 
than IO'^tq. (ii) Second, we quenched the system to 
T — 0.55e/kB below the melting and equilibrated it in 
a time interval longer than IO^tq. In this intermediate 
quenching, crystalline configurations were realizd^. (iii) 
Third, we further quenched the system to the final tem- 
perature T = 0.015e/fcB and annealed it in a time inter- 
val longer than 10'^ tq. There were almost no structural 
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changes in this step, (iv) Fourth, the periodic boundary 
condition along the x, y, z axises were changed such that 
the boundaries perpendicular to each axis became free 
surfaces. This change instantaneously caused a small de- 
compression and an expansion along each axis by about 
2%. We then annealed the system for a time longer than 
10"^To with T held fixed. At this low temperature, there 
was no detachment of the particles from the free surfaces 
due to the attractive interaction. After this procedure, 
we chose the particle positions in the boundary layers as 
the initial pinning points and introduced the spring 
potential in Eq.(5) between the bound particles and the 
pinning centers and the Nose-Hoover thermostats of the 
boundary layers. We set t = after these procedures. 



B. Bond breakage, disorder variable, and stress 

In dense particle systems, configuration changes can 
be visualized with the method of bond breakage^. For 
each particle configuration at a time a pair of particles 
z G a and j £ /3 is considered to be bonded if 



ry(i) = \r,{t)-r^[t)\ < A^a^p. 



(10) 



We set Ai = 1.3, for which Aiaap is slightly longer than 
the first peak distance of the pair correlation function 
dapir). For our fee crystals for c = 1, the corner- face dis- 
tance is equal to 1.088(T2 (which is slightly shorter than 
the minimum distance 2^/^(72 = 1.122(72 of the Lennard- 
Jones potential). The corner-corner distance or the lat- 
tice constant is then a = 1.088 x 2^/^a2 = 1.539ct2. Thus, 
for Ai = 1.3, the 12 nearest particle pairs at corner-face 
and face-face positions in a unit cell are bonded and the 
pairs at corner-corner positions are not bonded. After a 
time interval At, the bond is regarded to be broken if 



j{t + At) > A^Ua 



(11) 



This definition of bond breakage is insensitive to the value 
of A2 as long as A2 is larger than Ai and A2aap is shorter 
than the second peak position of gai^if)'^'^ ■ In our sim- 
ulation, the particle pairs with a broken bond in this 
definition underwent a relative motion mostly of order 
<Jai3- The broken bond number in time interval Ai] 
will be denoted by ANi,{t) as a function of time t at fixed 
At. The number of the particles with broken bonds in 
this time interval is twice of ANi,{t). 

Next, to examine the structural order, we introduce 
Steinhardt order parameter^^HlZl^ 
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which are defined for each unbound particle j. Here, 
Ygm{f) are the spherical harmonics functions of degree 
£ with —£< m < £, depending on the direction of a 
vector r. In our system, the local crystal structures are 
mostly fee and we choose £ = 6. In Eq.(12), we set 



Tjk (t) = Vj {t)—rk (t) , where the particle j G a at position 
rj{t) and the particle fc g /3 at position rk{t) are bonded 
{rjk < Aiaap)- The number of the particles bonded to 
the particle j is written as nl{t). This criterion of bonded 
particles is the same as in Eq.(lO), so nl{t) is mostly 
equal to 12 in fee crystal regions. We further introduce 
a disorder variable for each particle j bj'^Hl 
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= ^ E E \iLit)^<iL{t)\', (13) 

""hV*^ fcebond m=-6 

where the average over the bonded particles j is taken. 
For particles around defects Dj{t) is in the range 
[0.1, 0.7], while for a perfect crystal Dj{t) is nearly zero at 
low T {Dj < 0.005 in fee regions for c = 1). Thus D^{t) 
represents the deviation from hexagonal crystalline order 
for each particle j. It is convenient to introduce the av- 
erage of Dj{t) over the A^ub unbound particles (j G ub). 



(14) 
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which represents the overall degree of disorder. We have 
originally devised a disorder variable in two dimension^^. 

During stretching, we also calculated the average of 
the zz component of the stress'^^. It is expressed as sum- 
mations over the unbound particles in the form. 
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where {xj,yj,Zj) are the particle positions, {ij,yj,Zj) 
are the velocities, rjk = rj — rj, = {xjk, Vjk, zjk), are the 
relative position vectors with rjk = \rjk\, and Va^ijjk) 
are the pair potentials with v'^p{r) — dvap{r) / dr. Here 

V is the volume in the bulk interior region, but we set 

V = (Lq — 2c?o)io neglecting its weak time-dependence. 
We expect that {(Jzz){t) is approximately equal to the 
stress acting on the top and bottom after smoothing out 
of the contributions from the rapid spring motions^. To 
be precise, the overall stress distribution is somewhat 
complicated in our system due to the clamping condition 
at the top and bottom and the free-boundary condition 
on the side walls. See inward deformations of the side 
boundaries at strain ~ 0.06 in Fig. 5 below. 



III. SIMULATION RESULTS 

A. Stacking faults at t = 

We mostly realized intrinsic stacking faults with se- 
quence ABC AC ABC at low levels of disorder in un- 
strained states and under stretching. Furthermore, with 
addition of structural disorder in dilute mixtures, we also 
realized twin faults with sequences ABC ABC B AC B A 
during crystallization (in the preparation process), where 
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(a) cf=0.0594 




FIG. 1: (Color online) Relatively disordered particles with 
Dj > Do forming stacking faults in fee crystals at t — 0, 
where we set [0,02/ <Ji, Dq) equal to (a) (1,1.2,0.05), (b) 
(0.95,1.2,0.03), (c) (0.02,1.4,0.03), and (d) (0.95,1.4,0.05). 
Particles are written as spheres with diameter a\ or (T2 . Their 
colors represent Dj according to the color bar. Stacking 
faults visualized are intrinsic ones extending throughout the 
cell in (a), while they are twin ones stemming from a grain- 
boundary- like disordered region in the upper part in (b). 
Twin and intrinsic ones are mixed in (c). Successive intrinsic 
ones form hep regions in (d). 



the middle single-particle layer C forms a defect plane. 
We did not encounter extrinsic stacking faults with se- 
quence ABCABCBABC, where the middle layer B is 
inserted. 

The disordered structures in the initial particle config- 
urations at i = have been created during the crystUiza- 
tion in the preparation process. In Fig.l, using the initial 
data of unstrained states, we display relatively disordered 
particles with Dj > J?o in four cases: (a) c = 1 and 
(Ta/cTi = 1.2, (b) c = 0.95 and 02/(71 = 1.2, (c) c = 0.02 
and (T2/C1 — 1-4, and (d) c — 0.95 and o'2/cri = 1.4. 
The lower bound Dq is 0.05 in (a) and (d) and is 0.03 in 
(b) and (c) , since the particles composing intrinsic stack- 
ing faults have slightly larger Dj than those composing 
twin faults. Relatively ordered particles with Dj < Dq 
are transparent in these pictures. Then we can visualize 
stacking faults, because the particles on these planes have 
higher values of Dj than in the background fee region. 

In the single component case, we realized one or two 
intrinsic stacking faults extending throughout the cell in 
many cases, as in (a). We could also realize initial config- 
urations without stacking faults for c = 1. In the upper 
part of (a), two parallel stacking faults meet with sepa- 



(a) cf=0.0594 


(b) ef=0.0609 











(C) ef=0.1200 (d) ef=0.1215 



FIG. 2: (Color online) Particles twice undergoing partial rela- 
tive displacements on the (111) plane for c = 1 at et = 0.0594 
in (a) (top left), 0.0609 in (b) (top right), 0.1200 in (c) (bot- 
tom left), and 0.1215 in (d) (bottom right). Particle colors 
represent Dj{t) according to the color bar. In (a) a per- 
fect fee stacking is realized. Arrows indicate a partial rela- 
tive displacement between the top and next layers given by 
0.8(a/6)[121] here, which takes place just after this moment. 
In (b) an intrinsic stacking fault just created is shown after an 
elapse of 30to, where Dj of the top and next layers have in- 
creased to a higher value (in green). The third layer remains 
to have the lowest fee value of Dj (in blue). In (c), after a 
long elapse of time, the same stacking fault is going to disap- 
pear with a second partial displacement (a/6) [1 12] indicated 
by arrows. In (d), a perfect fee stacking is recovered. Data in 
Fig. 1(a) and Figs. 2-7 are from the same simulation run. 



ration of two particle layers. Partial dislocation cores are 
along the junction line. These stacking faults are parallel 
to the X axis and make an angle of O.SOtt with respect to 
the z axis. In (b), five twin faults with single-particle 
thickness can be seen. Four of them constitute two pairs 
of parallel twin faults stemming from a grain-boundary- 
like disordered region in the right upper part, while the 
bottom one is unpaired. For example, the two twin faults 
in the left upper part are paired as ABCBACABC with 
two C being visible. It is worth noting that paired twin 
faults can easily be produced from grain boundaries in 
nanocrystal^. 

In (c) and (d) of Fig.l, the size ratio 02! (J\ is increased 
from 1.2 to 1.4. See Fig. 6 of our previous paper^SJ for 
other snapshots of the particle configurations with the 
same parameter values, where disordered particles with 
Dj > 0.2 were shown to form tube-like or plate-like ag- 
gregates. At c = 0.02 in (c), twin and intrinsic stacking 
faults are mixed, where a small number of the large parti- 
cles strongly disturb the fee crystalline order of the small 
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FIG. 3: (Color online) Left: Potential energy U\{t) averaged 
over 50 particles on the A-th layer below the stacking fault in 
Fig. 2 emerging at t — ti ~ 20to with ti = 198Qro (top) and 
disappearing att — t2 ~ 20ro with t2 — 4000ro (bottom). This 
stacking fault is composed of the layers A = and —1. Right: 
Displacement Af\{t) averaged over the same 50 particles at 
the formation (top) and at the disappearance (bottom) of the 
stacking fault. For A = 0, the maximum velocity (~ 0.1(T2/ro) 
is much slower than the edge velocity (~ 4(j2/ro). 
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FIG. 4: (Color online) Average potential energy U\{t) and 
average disordered variable D\{t) (inset) for A = 0,3, and 6 
over a long time interval. The layer element A = undergoes 
two partial displacements at t ~ 2000ro and 4000to, where 
D\{t) increases or decreases steeply for A = 0. Another slip 
takes place at t ~ 2500ro apart from these layer elements, at 
which the elastic energy density decrease suddenly for A = 3 
and 6 but does not change for A = 0. 



particles. In (d), intrinsic stacking faults appear succes- 
sively to form hep regions. The local composition in the 
hep regions is about 0.95, while that of the (transparent) 
fee regions is 0.995. In this case, a small number of the 
small particles tend to be expelled from the fee and hep 
crystal region&SS^. 

Here we introduce the surface free energy for planar 
defects per unit aregP. In our simulation, it is well-defined 
in unstrained states with small overall disorder, though 
its meaning becomes not clear in highly strained states 
(as will be the case in Fig. 3). That is, at low T, it is 
given by the excess potential energy per unit area, since 
the entropic contribution is small. Thus we define the 
potential energy of the particle j £ a as 

U3 = \Y.''cp{\r,-rk\), (16) 

k 

where fc g /3 represent the surrounding particles. As well 
as Dj, the values of Uj supported by the particles in 
stacking faults are higher than the background fee value 
Uicc (— ~7.40e here). For intrinsic stacking faults in 
the panel (a) in Fig.l, the excess potential energy Uj — 
Uicc is about 0.15e on the two stacking- fault layers and 
is about O.le on the next layers. Since the areal particle 
density is about 0.92(t^^, we estimate the surface energy 



for intrinsic stacking faults as 

7i„t -0.5e/a2. (17) 

Note that the hep structure is a succession of intrinsic 
stacking faults as in the panel (d) of Fig.l. Therefore, 
the hep structure should have a higher energy than the 
fee structure roughly by O.le per particle for the Lennard- 
Jones potential. On the other hand, for twin faults in the 
panel (b), Uj — C/fcc is about 0.16e on the stacking- fault 
layer and is about 0.07e on the next layers. Thus the 
surface energy for twin faults is estimated as 

7twi„ - O.Se/ai (18) 

In our previous papei^^Slj we have estimated the surface 
energy from twin faults in its Fig. 3. 

If we set (T2 = 2 A and e = 300/cb in Eq.(17), we have 
7int ~ 50m J/m^ and 7twin ^ 30m J/m^. It is worth not- 
ing that these estimated values are of the same order as 
those for Cu and A^^!^ jj^ such nanocrystalline metals, 
partial dislocations and stacking faults have been widely 
observed. Full dislocations are more frequently observed 
for metals such as Al and Ni with larger stacking fault 
energies (> lOOmJ/m^JpHl. For hard -core colloids, the 
stacking fault energy is nearly zero, where coexistence of 
fee and hep structures have been observed^^^. 
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B. Elementary stacking- fault motions 

In fee crystals, an intrinsic stacking fault appears with 
a relative displacement bi of type (a/6) (121) taking place 
between two adjacent {111} plane Ji'. The lattice con- 
stant a is close to 2^/^(T2 for c = 1 here. A stacking fault 
disappears upon occurrence of a second relative displace- 
ment 62 of type (a/6)(211). The vector sum results in 
a full Burgers vector 6 = 5i + 62 of type (a/2) (110). 
In Figs. 2-7, we examine the stacking fault dynamics in 
various aspects using a single simulation run for c — 1. 

In Fig. 2, we illustrate an example of succession of these 
two-step processes for c = 1. The first and second par- 
tial displacements occur at i = 2000ro and 4000to with 
a duration time of order 20to. Here Dj on the stack- 
ing fault increases after the first slip and again decreases 
after the second slip. In this case, however, the size of 
the first relative displacement depicted in (a) is 80% of 
|5i| = a/^/2 = 0.72tT2 because of the elastic stress still 
remaining below the stacking fault, while the size of the 
second one in (c) is very close to \b2\ = a/^/2. 

To examine this slip dynamics in more detail, let us 
consider {111} layers close to the intrinsic stacking fault 
treated in Fig. 2. We pick up particles in the A-th layer be- 
low the stacking fault, whose positions {x, y, z) are speci- 
fied by < (z-f0.636y)/o-2 + A-14.5 < 1, 4 < x/a2 < 13, 
9 < y/(T2 < 13, and 3 < 2;/ct2 < 10. The number of these 
particles N\ is close to 50. The stacking fault itself con- 
sists of the layers with A = and —1. In Fig. 2, we have 
visualized the layers with A = —1,0, and 1. Here we take 
the average of the potential energy Uj in Eq.(16) over 
these particles to obtain 

Ux{t)= U,{t)/N^. (19) 

J G layer A 

The left panels of Fig. 3 display U\{t) versus t for A = 
0, 3, and 6 at the first and second slips. For A = 0, it 
exhibits a peak at each slip, where the peak height is 
higher by 0.2e than its stating value Uo{ti) or Uo{t2)- 
The resultant energy barrier for the (first and second) 
partial displacements per unit area is estimated as 

A7pcak - 0.2ea2^, (20) 

which is of the same order as 7int in Eq.(17)-. Right 
after the first displacement, U\{t) drops below tj\{ti) 
by 0.05e for A = and by 0.03e for A = 3, while it is 
nearly unchanged for A = 6. The elastic energy stored at 
t = ti \s not released for A = 6 in this case. In contrast, 
in unstrained states in Fig.l, the particles on stacking 
faults have higher Uj (by O.le here) than those in the 
fee regions, leading to a well-defined 7int. On the other 
hand, the local kinetic energy is enhanced only during 
the burst periods. That is, it increased by 0.05e for A = 
and by O.Ole for A = 6 per particle right after the 
two displacements, but the kinetic-energy peaks rapidly 
decayed on a time scale of 20ro (see K{t) in Fig. 5 and 
the remark (1) in the last section). 



The right panels of Fig. 3 display the time-dependent 
displacement of the A-th layer defined by 

ArA(t)= E |r,(t)-r,(<o)|/^A, (21) 

J tz layer A 

where to — ti — 1980ro or t2 — 4000to as in Fig. 3 and 
the average is taken over the N\ particles picked out in 
Eq.(19). In the first slip, the displacement at A = — 1 is 
opposite to that at A = with a nearly equal magnitude. 
In the second slip, the displacement at A = is over a 
distance of 0.38(72, while that at A = — 1 makes an angle 
of 37r/4 with respect to that at A = with a magnitude of 
0.45(72 • Then the relative displacement vector is nearly 
equal to 62 — (a/6) [211] as in the panel (c) of Fig. 2. 
The maximum of the sliding velocity of the layer A = 
is about 0.12(72/ro both for the first and second partial 
displacements, so the maximum of the relative velocity 
between the two layers is 0. 24(72 /tq. We shall see that 
the propagation velocity of a stacking fault edge is about 
4(72/ro in Fig. 7. 

In the left panels in Fig. 3, Uo{t) is lowest right after the 
first slip but is largest right before the second slip among 
A = 0, 3, and 6, owing to inhomogeneous accumulation 
of the elastic energy. In Fig. 4, we thus presents time- 
evolution of U\{t) on a long time interval in the same 
simulation run. In addition to the two slips in Figs. 2 
and 3, we notice occurrence of another slip at i ^ 2500to 
away from the layer elements under consideration. Re- 
markably, at this remote slip, the elastic energy density 
decreases suddenly by 0.05e for A = 3 and 6, but it is 
unchanged for A = 0, resulting in C7o(i) > ^^{t) > U(,{t) 
for t > 2500to. In the inset of Fig. 4, we also display the 
average disorder variable D\(t) defined by 

D,{t)^ D,{t)/N^. (22) 

jG layer A 

For A = O, D\{t) increases at t = 2000ro and decreases 
aX t = 4000to, while it remains to be the background fee 
value for A = 3 and 6. 



C. Stacking faults as slip elements 

In Fig. 5, we follow time-evolution of four averaged 
quantities during stretching in the strain range Q < it < 
0.2 for c = 1, where the initial state is in the panel (a) of 
Fig.l. In its top panel, we plot the average stress {<Jzz){t) 
in Eq.(15) and the broken bond numbers ANi){t) defined 
by Eqs.(lO) and (11) with At — IOtq. In the initial stage 
< it < 0.02, it approximately holds the linear elastic 
relation, 

{azz){t) = Eit. (23) 

where the coefficient E is about GOe/crf. Though the 
top and bottom boundaries are clamped here, E is close 
to Young's modulus 3GK/{K -t- G/3) ~ 3G, where the 
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FIG. 5: (Color online) Time evolution of averages over the 
unbound particles for c — 1 under stretching: {azz){t), U{t), 
K{t), and D{t) from above. Broken bond numbers ANt{t) are 
also shown (top), where At = lOro = 0.0003/e and A2 = 1.31. 
Drops of {cr^^)(i) and U{t) and bursts of ANi{t) and K{t) 
occur simultaneously at slip events. Bottom: Dj{t) in a plane 
(± the X axis) at eight points (a)- (h) marked on the curve of 
D{t), where stacking faults are all intrinsic ones. Initial one in 
(a) disappears in (b), while new ones emerge in (c), (d), and 
(e). In (f) one disappears at f ~ 4000to, which corresponds 
to the second partial displacement in Figs. 2 and 3. In (h) a 
rhcp region is realized with sequence ABC AC AB ABAC B. 



shear modulus G is about 20e/iT2 £^nd the bulk modulus is 
about lOOe/crf for c = 1. At each slip, {(Jzz){t) suddenly 
decreases and AiVf,(t) exhibits a sharp peak with heights 
mostly about 300 (^ (Lo/o'2)^)- The maximum in this 
run is about 700 at et = 1.2. These bursts last on a time 
scale of order IOtq. The resultant plastic strain is of order 
/S.Ni,[t) / N^\y and the typical stress drop is estimated as 



-r I'T' 10*11 



£f=0. 05994 
Disorder variable: ABC 




gf= Q.Q599 4 




FIG. 6: (Color online) Formation of an intrinsic stacking fault 
as a slip for c = 1. In the left top panel, particles are dis- 
played at et = 0.05994, whose bonds are broken subsequently 
in < £t - 0.05994 < eAt. We set At = 50ro, Ai = 1.3, 
and A2 = 1.4. Three fee principal axes are written (in white). 
Other panels show particles in the region 9.93 < x/a2 < 11.03 
on a plane (_L the x axis). Right top: Those after the slip 
around the two-particle layer (in red) with large opposite dis- 
placements, where the white bordered region in the right bot- 
tom panel is enlarged. Bottom: Those at it = 0.05994 (left) 
and at et — 0.06144 (right) before and after the stacking fault 
formation, respectively. Colors in the left top and bottom 
represent Dj{t) according to the color bar. 



which is consistent with the top panel of Fig. 5. 

The two middle panels of Fig. 5 give the potential and 
kinetic energies averaged over the unbound particles. 



jk^ ub 
jGub 



(25) 
(26) 



In defining the average potential, we should note that 
the particles on the free side boundaries have potential 
energies higher than those in the interior by 3 — 5e. Thus, 
in the summation of U (t) , we have removed the particles 
near the side boundaries with distances shorter than 3(72 . 
The number N^-^^ in Eq.(25) is then that of the unbound 
particles with larger separation from the side boundaries. 
After a slip at t = to, U (t) grows roughly as 



Act - GANb{t)/Nub - e/cr: 



(24) 



AU{t) = U{t) - U{to) - Ee^{t - tof/2. 



(27) 
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gf=0. 60003 £f=0. 06006 £f=0. 06003 ^/=0.06006 




FIG. 7: (Color online) Birth and growth of an intrinsic stacking fault for c — 1. Left four panels: Particles with their bonds 
being broken in successive time intervals neAt < et — 0.60003 < (n + l)eAt (n — 0,1,2,3) with At = to, Ai = 1.3, and 
A2 = 1.31. Colors represent Dj according to the color bar. Arrows (in yellow) indicate the direction of the fastest expansion 
of the partial dislocation curve. Right four panels: Particles on a plane in the range 12.4 < x/a2 < 13.5 with displacements 
multiplied by four. Red ones are those with their bonds broken in each time interval with Ai = 1.3 and A2 = 1.31. The edge 
velocity is 4(T2/ro. The partial dislocation core is located behind the edge within the ellipse in each panel. 



At each slip, a large fraction of the elastic energy is even- 
tually absorbed by the thermostats at the top and bot- 
tom very effectively because of short tnh and the small 
system size in our simulation. However, it is partly used 
to create the new stacking fault. Formation of a stacking 
fault with area Lq yields an increase in U {t) of order, 

AU,t ~ 7stio/A^ub ~ 0.02e. (28) 

On the other hand, the average kinetic energy K(t) ex- 
hibits sharp peaks smaller than 0.04e decaying to zero 
rapidly, so K{t) = 3kBT/2 = 0.0225e holds except for 
the burst periods. Thus heating is negligible. Further- 
more, D{t) in Fig. 5 increases or decrease in a stepwise 
manner upon appearance or disappearance of a stack- 
ing fault. This coincidence can be seen from side views 
of stacking faults in the eight bottom plates (a)-(h) in 
Fig. 5. Between (e) and (h) a stacking fault disappears, 
which corresponds to the second partial displacement in 
Figs. 2 and 3. In (h) we can see formation of a random 
sequence of stacking (a random hep region). 

In Fig. 6, we illustrate a slip event for c — 1, where an 
intrinsic (111) stacking fault is created in a time interval 
< it - 0.05994 < iAt = 0.0015. We set At = 50to, 
Ai — 1.3, and A2 ~ 1.4. The system length is in- 
creased only by 0.02(T2 in this time interval. In the left 
top, the displayed particles are just before the slip, so 



they still have the background fee value of Dj except 
for those near the free boundaries. Here the crystal ori- 
entation can be known from the crystal principal axes 
depicted, where [OlO] is parallel to the y axis, while [001] 
and and [100] make angles of ±7r/4 with respect to the 
z axis in the xz plane. Thus the slip plane makes an 
angle of cos^^(l/\/3) = O.SOtt with respect to the free 
surfaces (parallel to the xz plane). The other panels 
give cross-sectional particle configurations in the region 
9.93 < x/(j2 < 11.03. The right top panel displays the 
layer stacking and the displacements at the slip. The 
bottom panels show that Dj change on the stacking fault 
plane after the slip event. 

The formation of a stacking fault takes place in a very 
short time (< IOtq) in our small cell. Thus, in Fig. 7, 
we give sequences of snapshots of its growth for c = 1, 
where the time width is taken to be At = tq. In each 
panel in the left, particles with their bonds being broken 
are displayed, where we set Ai = 1.3 and A2 — 1.31. The 
direction of the fastest expansion of the stacking fault is 
marked by arrows. The expansion velocity of the edge is 
of order 4(72 /to, which is slightly smaller than the trans- 
verse sound velocity cj_ = (G/TO2n)^/^ ^ 5a2/To- Note 
that the actual particle velocity undergoing a slip is of 
order 0.1cr2/To), as already demonstrated in Fig3, and is 
much slower than the edge-expansion velocity. Thus the 
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values of Dj of the displayed particles (except those near 
the free boundaries) are still close to their background 
fee value. In the four right panels, we visualize cross- 
sectional particle configurations around the expanding 
stacking fault, where it is at its birth in the upper left 
panel, expands through the cell, and colloids with the 
free boundary in the bottom panels. Here, because A2 
is only slightly larger than Ai, the edge represents an 
inception of the bond breakage or a precursor of the par- 
tial dislocation. The relative displacement approaches its 
saturation length at the partial dislocation within 

the ellipses in the right panels in Fig. 7. The particles 
within these ellipses mostly have 11 bonds in our defini- 
tion in Eq.(lO). This precursor zone increases in time to 
have a width of 4-8 particles in the present case. In addi- 
tion, we note that the instantaneous particle velocities rj 
around the edge are rather random with large deviations 
because of their rattling motions in the crystal structure. 
Their time averages over tq given by J^'^'^° dt'rj{t')/To 
become proportional to the displacement vectors in the 
panels. 

We finally examine how stacking faults are removed 
by a collective partial displacement in applied strain. 
We start with the initial configuration (b) in Fig.l for 
c — 0.95, where the initial stacking faults were twin ones 
but an intrinsic stacking fault was created subsequently. 
Here the degree of disorder is considerably larger than in 
the case of c = 1. The left top panel of Fig. 8 displays 
particles whose bonds are just before breakage. They 
are mostly on the disappearing intrinsic stacking fault 
plane, but those in the projected part are in the grain 
boundary-like region in the upper part in the panel (b) 
of Fig.l. They are thus relatively disordered before the 
slip. The right top panel gives a cross-sectional profile 
of the stacking fault just before its disappearance, where 
the displacement vectors are from the present to future 
positions after a time elapse of At = SOtq. In the bot- 
tom plates, the values of the disorder variable are dis- 
played in the plane before and after the second partial 
displacement. The stacking CBCA (_L [111]) in the left 
disappears in the right by the registry change: B ^ A, 
C ^ B, and A — > C as marked in the left bottom corners. 



IV. SUMMARY AND REMARKS 

In Lennard- Jones particle systems, we have studied 
the dynamics of stacking faults with the aid of the 
bond breakage, the disorder variable Dj{t), and the 
particle displacements at T = O.OlSe/fcs. Hereafter, we 
summarize our main results together with comments, 
(i) We can select relatively disordered particles with Dj 
larger than a threshold Dq. For appropriate Dq, we can 
visualize stacking faults because the particles belonging 
to them have higher values of Dj than those in the 
fee crystal region. In Fig.l, we have shown the initial 
structural disorder for four cases of the composition c 
and the size ratio (72/(Ti to find intrinsic stacking faults. 




£f=0.0765 ef=0.0765 
Disorder variable: 




FIG. 8: (Color online) Removal of a preexisting intrinsic 
stacking fault by a subsequent partial displacement for c = 
0.95. Colors represent Dj according to the color bar. Top 
left; Particles whose bonds are just before breakage at the 
second displacement, so they mostly have relatively large Dj 
before the slip. Top right: Particles in a plane in the range 
3.9 < x/a2 < 5.0. Arrows represent displacements multiplied 
by 1.5 in a time interval with width At = 50t. Red particles 
have broken bonds with A2 — 1.35. Bottom: Dj of particles 
on the plane before the slip at it = 0.0765 (left) and after the 
slip at et = 0.0780 (right). 



twin faults, and hep regions. They indicate emergence 
of a vast variety of of the structural disorder sensitively 
depending on c and a2/<Ji- 

(ii) We have shown that plastic deformations in fee 
crystal can be achieved by partial displacements of 
close-packed planes, as illustrated in Fig. 2. We have 
examined rapid time-development of the potential 
energies and the positions of the particles on a stacking 
fault and those on layers nearby at the first and second 
partial displacements in Fig. 3. The potential energy 
of the particles increases by 0.2e when they are on an 
emerging or disappearing stacking fault. The relative 
particle velocity at stacking faults is of order 0.24cr2/To 
and is much slower than the transverse sound velocity 
c± ^ 5(72 /to- We have detected the position-dependent 
elastic energy stored during the elastic periods and its 
collective release upon slips in Fig. 4. 

(iii) We have shown time-evolution of some physical 
quantities over long times during stretching in Fig. 5. 
Upon appearance and disappearance of stacking faults, 
the stress exhibits strong intermittent fluctuations with 
collective bond breakage. The average elastic energy 
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accumulated in the cell is suddenly decreased upon slips, 
(iv) Apperance of an intrinsic stacking fault has been 
illustrated in Fig. 6. The formation of stacking faults 
occurs very rapidly, so Fig. 7 has given the time-evolution 
of a stacking fault on a time scale of tq. Remarkably, 
the edge of the stacking fault expands with a velocity 
close to the transverse sound velocity c±. As a result, 
there appears a precursor zone between the edge and the 
partial dislocation, whose width increases in time and is 
in a range of 7-10 particles in our case. Disappearnce of 
an intrinsic stacking fault has been illustrated in Fig. 8, 
where intrinsic stacking faults proliferate to form a rhcp 
sequence on a long time scale, however. 

We further make critical remarks as follows: 

(1) Because of the choice Ai = 1.3 in the definition of 
bonds in Eq.(lO), we have picked up the 12 nearest neigh- 
bors as bonded particles in fee crystal. The disorder vari- 
able Dj in Eq.(13) has then aquired the meaning of the 
deviation of hexagonal order in the fee crystal structure. 

(2) In our simulation, we have realized intrinsic stack- 
ing faults during stretching. However, in previous ex- 
periments on nanocrystal plasticitj^, paired twin faults 
like BACBABCABCACB have been observed to ap- 
pear from grain boundaries, where the second and last 
A form twin faults. 

(3) This paper has treated elementary dynamics of stack- 
ing faults in a small system with free boundaries. Thus 
stacking faults have appeared from the free side bound- 
aries. Larger system sizes are needed to examine large- 
scale plastic events. In particular, we should examine 
plastic deformations in polycrystals with many grains, 
where birth and growth of partial and full dislocations 
are strongly influenced by the grain boundaries. With de- 
creasing the grain size, the grain-boundary sliding should 
also come into play. 

(4) As indicated by Fig.l, the structural disorder becomes 
very complicated with increasing the composition of the 



second component. This disordering effect is much in- 
tensified for larger size ratios between the two particle 
species. Thus we should study the effects of increasing 
disorder in binary mixtures with size dispersity. 

(5) To understand the slip behavior in nanocrystals, 
Swygenhoven et al^ proposed to use a generalized stack- 
ing fault energy as a function of the particle displace- 
ment, because there are potential barriers for the forma- 
tion of partials. From the left panels of Fig. 3, we rec- 
ognize that the generalized stacking fault energy should 
furthermore depend on the local elastic strain. 

(6) We have attached Nose- Hoover thermostats at the top 
and bottom layers in the cell and have applied stretching 
by pulling the particles in the top layer. We should ex- 
amine this method in more detail by varying the param- 
eters such as the characteristic time tnh and the strain 
rate e. In this paper, tnh is very small {— 0.072ro) and 
the thermal equilibration after plastic events is very fast 

20to). However, for longer tnh, the thermostats be- 
come less efficient and the bulk thermal relaxation be- 
comes longer, resulting in heating in the cell. In larger 
systems, plastically deformed regions should emit sound 
waves, which propagate throughout the solid region and 
are reflected at the boundaries. This acoustic effect 
should be studied in future. 
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